Strategic Analysis of Micromobility Dock Deployment

Author
Affiliation

Joshua Sweren

Georgetown University

Code
#install.packages(c("osmdata", "sf", "dplyr", "ggplot2", "rnaturalearth"))
library(osmdata)
Code
library(sf)
Code
library(dplyr)
Code
library(ggplot2)
library(rnaturalearth)
library(osmextract)
Code
library(tigris)
Code
library(mapview)
library(lehdr)
library(viridis)
Code
library(viridisLite)
library(spatstat)
Code
library(units)
Code
library(spdep)
Code
library(terra)
Code
library(raster)
Code
library(tictoc)
Code
library(tidycensus)
library(leaflet)
Source: Article Notebook

Plotting the Docks

Code
dc_bb <- c(xmin = -77.45, ymin = 38.70, xmax = -76.80, ymax = 39.10)
docks_raw_dc <- opq(bbox = dc_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Capital Bikeshare") |>
  osmdata_sf()

docks_dc <- docks_raw_dc$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))
Code
hull_dc <- docks_dc |> 
    st_union() |> 
    st_convex_hull()
hull_buffered_dc <- st_buffer(hull_dc, dist = 2500)
hull_buffered_wgs84_dc <- st_transform(hull_buffered_dc, 4326)
docks_dc <- st_intersection(
    st_transform(docks_dc, 4326), 
    hull_buffered_wgs84_dc)
Code
mapview(hull_buffered_wgs84_dc, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_dc, col.regions = "red")
Code
nyc_bb <- c(xmin = -74.15, ymin = 40.55, xmax = -73.70, ymax = 40.95)
docks_raw_nyc <- opq(bbox = nyc_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Citi Bike") |>
  osmdata_sf()

docks_nyc <- docks_raw_nyc$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))
Code
hull_nyc <- docks_nyc |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_nyc <- st_buffer(hull_nyc, dist = 500)
hull_buffered_wgs84_nyc <- st_transform(hull_buffered_nyc, 4326)
docks_nyc <- st_intersection(
  st_transform(docks_nyc, 4326),
  hull_buffered_wgs84_nyc
)
Code
mapview(hull_buffered_wgs84_nyc, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_nyc, col.regions = "red")
Code
bos_bb <- c(-71.20, 42.20, -70.90, 42.45)
docks_raw_bos <- opq(bbox = bos_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Bluebikes") |>
  osmdata_sf()

docks_bos <- docks_raw_bos$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))
Code
hull_bos <- docks_bos |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_bos <- st_buffer(hull_bos, dist = 2500)
hull_buffered_wgs84_bos <- st_transform(hull_buffered_bos, 4326)
docks_bos <- st_intersection(
  st_transform(docks_bos, 4326),
  hull_buffered_wgs84_bos
)
Code
mapview(hull_buffered_wgs84_bos, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_bos, col.regions = "red")
Code
sf_bb <- c(-122.55, 37.70, -122.35, 37.85)
docks_raw_sf <- opq(bbox = sf_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Bay Wheels") |>
  osmdata_sf()

docks_sf <- docks_raw_sf$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))
Code
hull_sf <- docks_sf |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_sf <- st_buffer(hull_sf, dist = 500)
hull_buffered_wgs84_sf <- st_transform(hull_buffered_sf, 4326)
docks_sf <- st_intersection(
  st_transform(docks_sf, 4326),
  hull_buffered_wgs84_sf
)
Code
mapview(hull_buffered_wgs84_sf, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_sf, col.regions = "red")
Code
chi_bb <- c(-87.80, 41.75, -87.55, 42.00)
docks_raw_chi <- opq(bbox = chi_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Divvy") |>
  osmdata_sf()

docks_chi <- docks_raw_chi$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))
Code
hull_chi <- docks_chi |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_chi <- st_buffer(hull_chi, dist = 2500)
hull_buffered_wgs84_chi <- st_transform(hull_buffered_chi, 4326)
docks_chi <- st_intersection(
  st_transform(docks_chi, 4326),
  hull_buffered_wgs84_chi
)
Code
mapview(hull_buffered_wgs84_chi, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_chi, col.regions = "red")

Plotting Other Data

Code
#census_api_key("5e8593d6e289e6ae04c33a019f4951ca2bdd9523", install = TRUE)
dc_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("DC", "MD", "VA"),
  geometry = TRUE
)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |======================================================================| 100%

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |======================================================================| 100%

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |======================================================================| 100%
Code
dc_bg <- st_transform(dc_bg, st_crs(hull_buffered_wgs84_dc))

dc_bg <- dc_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

dc_bg_clipped <- st_intersection(dc_bg, hull_buffered_wgs84_dc)
Code
dc_bg_clipped <- dc_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

dc_bg_density <- dc_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )
Code
max_density <- 20000
dc_bg_density$pop_density_plot <- pmin(dc_bg_density$pop_density, max_density)
Code
dc_bbox <- st_bbox(hull_buffered_wgs84_dc)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_dc <- opq(dc_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_dc))
Code
pal <- leaflet::colorNumeric(viridis(20), domain = dc_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = dc_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_dc,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_dc,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )
Code
nyc_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("NY", "NJ"),
  geometry = TRUE
)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |======================================================================| 100%

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |======================================================================| 100%
Code
nyc_bg <- st_transform(nyc_bg, st_crs(hull_buffered_wgs84_nyc))

nyc_bg <- nyc_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

nyc_bg_clipped <- st_intersection(nyc_bg, hull_buffered_wgs84_nyc)
Code
nyc_bg_clipped <- nyc_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

nyc_bg_density <- nyc_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )
Code
max_density <- 75000
nyc_bg_density$pop_density_plot <- pmin(nyc_bg_density$pop_density, max_density)
Code
nyc_bbox <- st_bbox(hull_buffered_wgs84_nyc)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_nyc <- opq(nyc_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_nyc))
Code
pal <- leaflet::colorNumeric(viridis(20), domain = nyc_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = nyc_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_nyc,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_nyc,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )
Code
bos_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("MA"),
  geometry = TRUE
)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |======================================================================| 100%
Code
bos_bg <- st_transform(bos_bg, st_crs(hull_buffered_wgs84_bos))

bos_bg <- bos_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

bos_bg_clipped <- st_intersection(bos_bg, hull_buffered_wgs84_bos)
Code
bos_bg_clipped <- bos_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

bos_bg_density <- bos_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )
Code
max_density <- 25000
bos_bg_density$pop_density_plot <- pmin(bos_bg_density$pop_density, max_density)
Code
bos_bbox <- st_bbox(hull_buffered_wgs84_bos)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_bos <- opq(bos_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_bos))
Code
pal <- leaflet::colorNumeric(viridis(20), domain = bos_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = bos_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_bos,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_bos,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )
Code
sf_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("CA"),
  geometry = TRUE
)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |======================================================================| 100%
Code
sf_bg <- st_transform(sf_bg, st_crs(hull_buffered_wgs84_sf))

sf_bg <- sf_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

sf_bg_clipped <- st_intersection(sf_bg, hull_buffered_wgs84_sf)
Code
sf_bg_clipped <- sf_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

sf_bg_density <- sf_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )
Code
max_density <- 25000
sf_bg_density$pop_density_plot <- pmin(sf_bg_density$pop_density, max_density)
Code
sf_bbox <- st_bbox(hull_buffered_wgs84_sf)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_sf <- opq(sf_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_sf))
Code
pal <- leaflet::colorNumeric(viridis(20), domain = sf_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = sf_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_sf,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_sf,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )
Code
chi_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("IL"),
  geometry = TRUE
)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%
Code
chi_bg <- st_transform(chi_bg, st_crs(hull_buffered_wgs84_chi))

chi_bg <- chi_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

chi_bg_clipped <- st_intersection(chi_bg, hull_buffered_wgs84_chi)
Code
chi_bg_clipped <- chi_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

chi_bg_density <- chi_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )
Code
max_density <- 30000
chi_bg_density$pop_density_plot <- pmin(chi_bg_density$pop_density, max_density)
Code
chi_bbox <- st_bbox(hull_buffered_wgs84_chi)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_chi <- opq(chi_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_chi))
Code
pal <- leaflet::colorNumeric(viridis(20), domain = chi_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = chi_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_chi,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_chi,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

Univariate Testing

Code
coords <- st_coordinates(docks_dc)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)
Code
lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 35.752, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.8444759043     -0.0012936611      0.0005596277 
Code
dock_values <- inv_dist

localI <- localmoran(dock_values, lw, zero.policy = TRUE)
docks_dc$localI <- localI[,1]

leaflet(docks_dc) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    radius = 4,
    color = ~viridis::viridis(20)[cut(localI, breaks = 20)],
    fillOpacity = 0.8,
    popup = ~paste0("Local Moran's I: ", round(localI, 3))
  )
Code
coords <- st_coordinates(docks_nyc)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)
Code
lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 48.697, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.6827748494     -0.0004344049      0.0001968368 
Code
coords <- st_coordinates(docks_bos)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 23.091, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.842640558      -0.003134796       0.001341649 
Code
coords <- st_coordinates(docks_sf)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 12.936, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.808442476      -0.009433962       0.003997522 
Code
coords <- st_coordinates(docks_chi)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 43.616, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.9584339511     -0.0011013216      0.0004839866 

Results

Metropolitan Area Moran’s I Statistic
Washington, DC \(0.8444759043\)
New York, NY \(0.6827748494\)
Boston, MA \(0.842640558\)
San Francisco, CA \(0.808442476\)
Chicago, IL \(0.9584339511\)

Hypothesis Test: Population Density

Code
proj_crs <- "EPSG:32618"
docks_proj <- st_transform(docks_dc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_dc, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
dc_bg_density_proj <- st_transform(dc_bg_density, 26918)
rtemplate <- rast(
  extent = ext(dc_bg_density_proj),
  resolution = 250,
  crs = crs(dc_bg_density)
)

pop_rast <- rasterize(
  vect(dc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)
Code
dc_bg_density_proj <- st_transform(dc_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(dc_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(dc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))
Code
nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  774 points
Average intensity 5.06e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 822.6608 x 638.4921 units

Original dummy parameters: =
Planar point pattern:  2974 points
Average intensity 1.94e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [47800, 525000]  total: 1.53e+09
Weights on data points:
    range: [47800, 263000]  total: 1.32e+08
Weights on dummy points:
    range: [47800, 525000]  total: 1.4e+09
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.512591e+01  1.913035e-04 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.512591e+01 4.795913e-02 -1.521990e+01 -1.503191e+01   ***
pop_im       1.913035e-04 5.669327e-06  1.801919e-04  2.024152e-04   ***
                  Zval
(Intercept) -315.39161
pop_im        33.74361

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.512591e+01  1.913035e-04 

Fitted exp(theta):
 (Intercept)       pop_im 
2.697132e-07 1.000191e+00 
Code
anova(m0, m1, test="LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   683.18 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
docks_proj <- st_transform(docks_nyc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_nyc, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
nyc_bg_density_proj <- st_transform(nyc_bg_density, 26918)
rtemplate <- rast(
  extent = ext(nyc_bg_density_proj),
  resolution = 250,
  crs = crs(nyc_bg_density)
)

pop_rast <- rasterize(
  vect(nyc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)
Code
nyc_bg_density_proj <- st_transform(nyc_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(nyc_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(nyc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))
Code
nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  2303 points
Average intensity 5.99e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568

Dummy quadrature points:
     100 x 100 grid of dummy points, plus 4 corner points
     dummy spacing: 215.4011 x 314.1755 units

Original dummy parameters: =
Planar point pattern:  5780 points
Average intensity 1.5e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568
Quadrature weights:
     (counting weights based on 100 x 100 array of rectangular tiles)
All weights:
    range: [8530, 67700]    total: 3.84e+08
Weights on data points:
    range: [11300, 33800]   total: 69200000
Weights on dummy points:
    range: [8530, 67700]    total: 3.15e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.237567e+01  2.121081e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.237567e+01 3.074223e-02 -1.243592e+01 -1.231541e+01   ***
pop_im       2.121081e-05 9.451471e-07  1.935836e-05  2.306327e-05   ***
                  Zval
(Intercept) -402.56236
pop_im        22.44181

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.237567e+01  2.121081e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
4.220043e-06 1.000021e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 4.3% (345 out of 8083) of the quadrature points 
Code
anova(m0, m1, test="LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1  346                          
2  347  1   436.56 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
docks_proj <- st_transform(docks_bos, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_bos, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
bos_bg_density_proj <- st_transform(bos_bg_density, 26918)
rtemplate <- rast(
  extent = ext(bos_bg_density_proj),
  resolution = 250,
  crs = crs(bos_bg_density)
)

pop_rast <- rasterize(
  vect(bos_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)
Code
bos_bg_density_proj <- st_transform(bos_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(bos_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(bos_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))
Code
nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  320 points
Average intensity 1.17e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 305.5909 x 310.2755 units

Original dummy parameters: =
Planar point pattern:  2940 points
Average intensity 1.08e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [12100, 94800]   total: 2.72e+08
Weights on data points:
    range: [19000, 47400]   total: 13400000
Weights on dummy points:
    range: [12100, 94800]   total: 2.59e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.418988e+01  9.188116e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.418988e+01 8.618496e-02 -1.435880e+01 -1.402096e+01   ***
pop_im       9.188116e-05 8.049611e-06  7.610421e-05  1.076581e-04   ***
                  Zval
(Intercept) -164.64449
pop_im        11.41436

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.418988e+01  9.188116e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
6.877236e-07 1.000092e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 4.4% (145 out of 3260) of the quadrature points 
Code
anova(m0, m1, test="LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1  146                          
2  147  1    105.4 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
docks_proj <- st_transform(docks_sf, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_sf, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
sf_bg_density_proj <- st_transform(sf_bg_density, 26918)
rtemplate <- rast(
  extent = ext(sf_bg_density_proj),
  resolution = 250,
  crs = crs(sf_bg_density)
)

pop_rast <- rasterize(
  vect(sf_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)
Code
sf_bg_density_proj <- st_transform(sf_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(sf_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(sf_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))
Code
nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  107 points
Average intensity 8.44e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743

Dummy quadrature points:
     32 x 32 grid of dummy points, plus 4 corner points
     dummy spacing: 398.6460 x 417.8386 units

Original dummy parameters: =
Planar point pattern:  805 points
Average intensity 6.35e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743
Quadrature weights:
     (counting weights based on 32 x 32 array of rectangular tiles)
All weights:
    range: [7200, 167000]   total: 1.27e+08
Weights on data points:
    range: [41600, 83300]   total: 8120000
Weights on dummy points:
    range: [7200, 167000]   total: 1.19e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.432761e+01  3.904106e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.432761e+01 1.743289e-01 -1.466929e+01 -1.398594e+01   ***
pop_im       3.904106e-05 1.500172e-05  9.638223e-06  6.844389e-05    **
                  Zval
(Intercept) -82.187281
pop_im        2.602438

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.432761e+01  3.904106e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
5.992334e-07 1.000039e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 0.66% (6 out of 912) of the quadrature points 
Code
anova(m0, m1, test="LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance Pr(>Chi)  
1    7                       
2    8  1   6.4363  0.01118 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
docks_proj <- st_transform(docks_chi, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_chi, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
chi_bg_density_proj <- st_transform(chi_bg_density, 26918)
rtemplate <- rast(
  extent = ext(chi_bg_density_proj),
  resolution = 250,
  crs = crs(chi_bg_density)
)

pop_rast <- rasterize(
  vect(chi_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)
Code
chi_bg_density_proj <- st_transform(chi_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(chi_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(chi_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))
Code
nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  909 points
Average intensity 9.85e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697

Dummy quadrature points:
     70 x 70 grid of dummy points, plus 4 corner points
     dummy spacing: 362.3008 x 745.5848 units

Original dummy parameters: =
Planar point pattern:  3482 points
Average intensity 3.77e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697
Quadrature weights:
     (counting weights based on 70 x 70 array of rectangular tiles)
All weights:
    range: [24600, 270000]  total: 9.21e+08
Weights on data points:
    range: [24600, 135000]  total: 1.08e+08
Weights on dummy points:
    range: [24600, 270000]  total: 8.14e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.419576e+01  9.193927e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.419576e+01 4.680314e-02 -1.428750e+01 -1.410403e+01   ***
pop_im       9.193927e-05 4.699411e-06  8.272859e-05  1.011499e-04   ***
                 Zval
(Intercept) -303.3079
pop_im        19.5640

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.419576e+01  9.193927e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
6.836891e-07 1.000092e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 10% (444 out of 4391) of the quadrature points 
Code
anova(m0, m1, test="LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1  445                          
2  446  1   269.89 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results

City Docks Dock Intensity (points/unit squared) pop_im coefficient Z-value Deviance diff
Washington, DC \(774\) \(5.06 \times 10^{-7}\) \(1.913 \times 10^{-4}\) \(33.74\) \(683.18\)
New York, NY \(2,303\) \(1.5 \times 10^{-5}\) \(2.121 \times 10^-5\) \(22.44\) \(436.56\)
Boston, MA \(320\) \(1.08 \times 10^{-5}\) \(9.189 \times 10^{-5}\) \(11.41\) \(105.4\)
San Francisco, CA \(107\) \(6.35 \times 10^{-6}\) \(3.904 \times 10^{-5}\) \(2.60\) \(6.436\)
Chicago, IL \(909\) \(3.77 \times 10^{-6}\) \(9.194 \times 10^{-5}\) \(19.56\) \(269.89\)

Hypothesis Test: Other Transit options

Code
stations_proj <- st_transform(stations_dc, 32618)
docks_proj <- st_transform(docks_dc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_dc, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)
Code
dist_im <- distmap(station_ppp)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   1112.3 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  774 points
Average intensity 5.06e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 822.6608 x 638.4921 units

Original dummy parameters: =
Planar point pattern:  2974 points
Average intensity 1.94e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [47800, 525000]  total: 1.53e+09
Weights on data points:
    range: [47800, 263000]  total: 1.32e+08
Weights on dummy points:
    range: [47800, 525000]  total: 1.4e+09
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-12.611735447  -0.001085935 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -12.611735447 0.0604945559 -12.730302598 -12.493168296   ***
dist_im      -0.001085935 0.0000453318  -0.001174783  -0.000997086   ***
                  Zval
(Intercept) -208.47720
dist_im      -23.95525

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-12.611735447  -0.001085935 

Fitted exp(theta):
 (Intercept)      dist_im 
3.332674e-06 9.989147e-01 
Code
stations_proj <- st_transform(stations_nyc, 32618)
docks_proj <- st_transform(docks_nyc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_nyc, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)
Code
dist_im <- distmap(station_ppp)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   1174.6 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  2303 points
Average intensity 5.99e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568

Dummy quadrature points:
     100 x 100 grid of dummy points, plus 4 corner points
     dummy spacing: 215.4011 x 314.1755 units

Original dummy parameters: =
Planar point pattern:  5780 points
Average intensity 1.5e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568
Quadrature weights:
     (counting weights based on 100 x 100 array of rectangular tiles)
All weights:
    range: [8530, 67700]    total: 3.84e+08
Weights on data points:
    range: [11300, 33800]   total: 69200000
Weights on dummy points:
    range: [8530, 67700]    total: 3.15e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-11.074435702  -0.001612625 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -11.074435702 3.398891e-02 -11.141052740 -11.007818663   ***
dist_im      -0.001612625 6.225954e-05  -0.001734652  -0.001490599   ***
                  Zval
(Intercept) -325.82498
dist_im      -25.90166

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-11.074435702  -0.001612625 

Fitted exp(theta):
 (Intercept)      dist_im 
1.550364e-05 9.983887e-01 
Code
stations_proj <- st_transform(stations_bos, 32618)
docks_proj <- st_transform(docks_bos, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_bos, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)
Code
dist_im <- distmap(station_ppp)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   263.95 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  320 points
Average intensity 1.17e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 305.5909 x 310.2755 units

Original dummy parameters: =
Planar point pattern:  2940 points
Average intensity 1.08e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [12100, 94800]   total: 2.72e+08
Weights on data points:
    range: [19000, 47400]   total: 13400000
Weights on dummy points:
    range: [12100, 94800]   total: 2.59e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-12.306690343  -0.001460137 

                 Estimate         S.E.     CI95.lo       CI95.hi Ztest
(Intercept) -12.306690343 0.0944239275 -12.4917578 -12.121622846   ***
dist_im      -0.001460137 0.0001187592  -0.0016929  -0.001227373   ***
                  Zval
(Intercept) -130.33445
dist_im      -12.29494

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-12.306690343  -0.001460137 

Fitted exp(theta):
 (Intercept)      dist_im 
4.521393e-06 9.985409e-01 
Code
stations_proj <- st_transform(stations_sf, 32618)
docks_proj <- st_transform(docks_sf, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_sf, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)
Code
dist_im <- distmap(station_ppp)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance Pr(>Chi)   
1    1                        
2    2  1    9.655 0.001888 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  107 points
Average intensity 8.44e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743

Dummy quadrature points:
     32 x 32 grid of dummy points, plus 4 corner points
     dummy spacing: 398.6460 x 417.8386 units

Original dummy parameters: =
Planar point pattern:  805 points
Average intensity 6.35e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743
Quadrature weights:
     (counting weights based on 32 x 32 array of rectangular tiles)
All weights:
    range: [7200, 167000]   total: 1.27e+08
Weights on data points:
    range: [41600, 83300]   total: 8120000
Weights on dummy points:
    range: [7200, 167000]   total: 1.19e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-1.362773e+01 -2.971307e-04 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.362773e+01 0.1466369755 -13.915136434 -1.334033e+01   ***
dist_im     -2.971307e-04 0.0001051373  -0.000503196 -9.106536e-05    **
                  Zval
(Intercept) -92.935177
dist_im      -2.826121

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-1.362773e+01 -2.971307e-04 

Fitted exp(theta):
 (Intercept)      dist_im 
1.206565e-06 9.997029e-01 
Code
stations_proj <- st_transform(stations_chi, 32618)
docks_proj <- st_transform(docks_chi, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_chi, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)
Code
dist_im <- distmap(station_ppp)
Code
m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")
Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   386.02 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
summary(m1)
Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  909 points
Average intensity 9.85e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697

Dummy quadrature points:
     70 x 70 grid of dummy points, plus 4 corner points
     dummy spacing: 362.3008 x 745.5848 units

Original dummy parameters: =
Planar point pattern:  3482 points
Average intensity 3.77e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697
Quadrature weights:
     (counting weights based on 70 x 70 array of rectangular tiles)
All weights:
    range: [24600, 270000]  total: 9.21e+08
Weights on data points:
    range: [24600, 135000]  total: 1.08e+08
Weights on dummy points:
    range: [24600, 270000]  total: 8.14e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-1.285494e+01 -8.244326e-04 

                 Estimate         S.E.       CI95.lo      CI95.hi Ztest
(Intercept) -1.285494e+01 5.571363e-02 -1.296414e+01 -12.74574347   ***
dist_im     -8.244326e-04 4.844097e-05 -9.193751e-04  -0.00072949   ***
                  Zval
(Intercept) -230.73239
dist_im      -17.01933

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-1.285494e+01 -8.244326e-04 

Fitted exp(theta):
 (Intercept)      dist_im 
2.613187e-06 9.991759e-01 

Results

City Docks Dock Intensity (points/unit squared) dist_im coefficient Z-value Deviance diff
Washington, DC \(774\) \(5.06 \times 10^{-7}\) \(-1.09 \times 10^{-3}\) \(-23.96\) \(1112.3\)
New York, NY \(2,303\) \(1.5 \times 10^{-5}\) \(-1.61 \times 10^{-3}\) \(-25.90\) \(1174.6\)
Boston, MA \(320\) \(1.08 \times 10^{-5}\) \(-1.46 \times 10^{-3}\) \(-12.29\) \(263.95\)
San Francisco, CA \(107\) \(8.44 \times 10^{-7}\) \(-2.97 \times 10^{-04}\) \(-2.83\) \(9.655\)
Chicago, IL \(909\) \(3.77 \times 10^{-6}\) \(-8.24 \times 10^{-4}\) \(-17.02\) \(386.02\)